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Abstract 

The long-time behavior of the current auto-correlation functions for the velocity, the shear stress 
and the heat flux is investigated in freely cooling granular gases. It is found that the correlation 
functions for the velocity and the shear stress have the long-time tails obeying r^'^/^, while the 
correlation function for heat flux decays as r"*-"^"^^^/^ exp(— C*r) with the dimensionless cooling rate 
C*, the spatial dimension d and the scaled time r in terms of the collision frequency. The result 
of our numerical simulation of the freely cooling granular gases is consistent with the theoretical 
prediction. 
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I. INTRODUCTION 



One of the central concerns in granular physics is to know the rheological properties of 
granular fluids [1]. It is believed that rapid granular flows for relatively dilute granular 

BLses can be described by a set of hydrodynamic equations derived from the kinetic theory 
, y, Q|. Most of theories assume the molecular chaos ansatz which can be only justified 
for dilute^gases. For granular gases with finite density, inelastic Enskog equation has been 
and the theoretical predictions well recover the results of simulations and 



used 

experiments [9|, 




It is, however, well known that correlated collisions cause significant 



changes in the form of constitutive relations for molecular fluids. It is important to include 
effects of correlated collisions for the description of granular fluids, though so far there are 
few such studies. 

In a recent paper, Saitoh and Hayakawa [ll|] demonstrate the relevancy of the kinetic 
approach to describe sheared granular fluids [5!]. In spite of their success, there are some 
unclear points. For example, they rely on the kinetic theory for uniform cooling granular 
gases, but this approach may not be valid for the description of shear flows. In fact, we 
have recognized the significant differences between freely cooling granular gases and sheared 
granular flows: (i) there is a long-range velocity correlation in freely cooling states 12|, but 
we do not find any evidence of the equal-time long-range correlation in sheared granular 
flows [12 1, (ii) the basic solution of the inelastic Boltzmann equation in sheared granular 



flows 
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151] is completely different from that of homogeneous cooling states 16|. 



In general, the long-time tails obeying t~'^^'^ with the time t and the spatial dimension 
d in the auto-correlation functions play important roles for elastic gases 17, Q, [3]. In 
fact, it is known that the transport coefficients in two-dimensional systems diverge in the 
thermodynamic limit 



20 



2l[ |. and they have the logarithmic singularity in the virial ex- 



pansions even in three dimensional cases [22] because of the long-time tails. Thus, if there 
are the long-time tails in granular fluids, we may need to change the transport coefficients 



In a 



determined by the inelastic Boltzmann equation jsj and Enskog equation 
recent paper, Kumaran 23| indicates the s upp resses of the long-time tail in sheared granu- 
lar flows as t"'^'^/^, while Ahmad and Puri [24;] suggest the existence of the long-time tail as 
with the scaled time r by the collision frequency in freely cooling granular gases from 
their two-dimensional simulation. Therefore, we need to clarify what is the true story of the 
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long-time tails in granular fluids and whether the logarithmic divergences of the transport 
coefficients exist in two-dimensional granular gases. 

In this paper, we discuss whether there are the long-time tails in freely cooling granular 
fluids. First, we demonstrate that the conventional long-time tails for the velocity auto- 
correlation function and the shear stress auto-correlation function exist, while we predict that 
the auto-correlation function for the heat flux decays exponentially from our consideration 
of the hydrodynamic fluctuations around the uniform cooling state. Second, we show the 
results of our extensive numerical simulations for two-dimensional freely cooling granular 
gases, which are almost consistent with our theory. Finally, we will discuss the non-trivial 
relations between the diffusion coefficients and the auto-correlation functions in granular 
ffuids. 



II. THEORETICAL ANALYSIS 

In this section, we calculate the behavior of the time correlation functions in the freely 
cooling system based on the method developed in ref. |18|. This section consists of four 
subsections associated with five appendices. 

A. Model 

The system considered in this paper consists of identical smooth and hard spherical 
particles with the mass m and the diameter a in the volume V . The position and the 
velocity of the i-th particle at time t are ri{t) and respectively. The particles collide 

instantaneously with each other with a coefficient of restitution e less than unity, where 
we do not consider the effects of oblique impacts and the collision speed dependence in e 



25l . |26| . |27| . |28| . |29| . [SOj. When the particle i with the velocity vi collides with the particle j 



with Vj, the post-collisional velocities v'i and v' j are respectively given by 

v'i = Vi-]^{1 + e){e- Vij)e, (1) 
1 

where e is the unit vector parallel to the relative position of the two colliding particles at 
contact, and vij = Vi — Vj. 
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B. Basic formulae 



Here, we introduce the correlation functions as 

1 ^ 

CD{to,t) = lim —y^(vi^{to)vi^{t + to))to, (3) 

V^oo V ■^"^ 
i=l 

C^{to,t) = lim i(J (to)J (t + to)),^, (4) 
Cx{to,t) = lim hj^{to)Jx{t + to))to, (5) 

V^oo V 

where Jr^(t) and J\(t) are the shear stress and the heat flux at time t, respectively. Here, 
to is the time when we start the measurement and (■ ■ ■ represents the ensemble average 
of possible conflguration at time to- In general, the currents Jriit) and Jx(t) consist of 
the kinetic part and the potential part, respectively. In this paper, we only consider the 
contribution from the kinetic part of the currents. This treatment is correct in dilute gases. 
For higher density cases, we need a more sophisticated method to include the contribution 
from the potential part, but the corrections only appear in the prefactor of coefficients at 
least for elastic gases 31, 3- To be consistent with the assumption that the kinetic part is 
dominant, our calculation is based on the hydrodynamic equations for dilute granular gases. 
Thus, the currents in Eqs. (jlj) and ([5]) are respectively approximated by their kinetic parts 
J^(t) and Jf (t) as 

J^{t) c:^ J^{t) = J2^^Ut)viy{t), (6) 

i 

Jxit) ^ J^{t) H'(^) - + 2)2^W) (7) 

i 

where T{t) is the temperature at time t. In order to simplify the calculation, we express 
these correlation functions as 

C„(to,t) = lim i(Arj,(i;i(to))Ja(t + to))to, (8) 

V^oo v 

where a = D,ri, X and 

joivi) = (9) 



j^{vi) = mvi^viy, (10) 
JxM = ^{mvf-{d + 2)T)v,,, (11) 
Jd = Jd{v^). (12) 



In the following argument, the calculation of the correlation functions will be separated 
into two steps. At the first step, we decompose the average in Eq. ([8]) into the partial aver- 
ages over spatially nonuniform systems where the particle 1 has a given velocity i'i(to) = ''^o 
and its initial position is constrained to the neighborhood by the smeared out proba- 
bility density W{ri{to) — vq). Thus, the decomposed initial weight function is given by 
6{vo — Vi{to))W{rQ — r'i(to)), where W{r) satisfies the normalization 



drW{r) = 1. (13) 

At the second step, we further average quantities over Vq and Vq in the given ensemble of 
granular fluids. 

Following this procedure, we rewrite Eq. as 

Ca{to,t) = ^i^^ j dro j dvoja{vo) 

X j dr{NJa{r,t + to)5ivo-Vi{to))Wiro-ri{to)))to, (14) 
where we introduce the microscopic current densities as 

JDir,t + to) = joMt + to))6{r - ri{t + to)), 

mvix{t + tQ)viy(t + to)S{r — ri{t + to)), 

i 

Jx{r, t + to) = J2^ {mv^it + to) -{d + 2)T{t + to)) ViJ{r -ri{t + to)). (15) 

i 

Here, we define the conditional average {F{r, t + to))^^ ^ of an arbitrary function F{r, t + to) 

as 

{NF{r, t + to)5{vo - v,{to))W{ro - ri(to))),„ = (5(^o - v,{to))W{ro - r'i(to))),„ {F{r, t + to)),„,e 

= h{to,Vo){F{r,t + to))^^^^^, (16) 

where fo(to,vo) is the one-particle velocity distribution function at the initial time to- The 
choice of foito^vo) is not trivial. If we are interested in the relaxation process starting from 
an equilibrium state, /o(io)'^o) can be the Maxwell-Boltzmann distribution. This choice, 
however, cannot be used for the case starting from the middle of a cooling state. If we 
assume that the base state is in a homogeneous cooling state (HCS), it is natural to adopt 
/o(^o? ^o) is the approximate solution of HCS as 

d/2 



fo{to,Vo) = riH 



m 



[27rTo(to) 



e-^ {l + a^S2{c')}, (17) 
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where To(to) is the initial temperature at to- In HCS, is approximately described as 

16(l-e)(l-2e) 
2 9 + 24rf + 8rfe-41e + 30(l -e)e2 



(18) 



determined from the first Sonine expansion with c = f A/m/2To(to) and S2{x) = — {d + 
2)x/2 + d{d + 2)/8 [16]. It should be noted that a more sophisticated method developed 



by Brilliantov and Poschel 33| which includes the higher expansion term af^ S^lc^) with 
5'3(x) = -a;V6 + ((i + 4)/4a;2-(d + 2)(d + 4)a;/8 + rf(rf + 2)(rf + 4)/48 does not give significant 



differences in the statistical properties of HCS from that by van Noije and Ernst 16|, 24 1. 
Hence, we neglect the higher order terms such as af^ for later discussions. 
With the aid of this conditional average, Eq. (|T4l) is reduced to 



Ca{to,t) = j dro J dvoja{vo)fo{to,Vo) j dr {Jo,{r,t + to))^^^^^. (19) 

We rewrite (^Ja{r,t + to))^^ ^ furthermore. For a = r],X, from Eq. (fT5|) . {Ja{r,t + to))to,c can 
be approximated by 

dvUv)fi{r,v,t). (20) 

Here, we replace ^i^iit + to) ~ 'r)5{vi{t + to) — v))^^ ^ by the local scaling distribution 
function |16[ 

-1 d/2 



fi{r,v,t + to) =n{r,t + to) 



_27rT(r,t + to 



e-^'{l + 02^2(5') + ^ap^p(c2)}, (21) 



p=3 



where c = ^y^m\v — u{r, t + to)p/ (2T(r, t + to)) It should be noted that 02 is different from 
and with I > 3 may be relevant for inhomogeneous cooling state. The replacement 
of ^{.fiit + to) — r)5{vi{t + to) — 'i'))(„ c flight be crucial. However, (i) in general, 
the velocity distribution function (VDF) can be expanded around the local Maxwellian, (ii) 
the transport coefficients and the velocity correlations can be determined from the lower 
moments of VDF and (iii) the long-time tails do not depend on the choice of fi. From Eq. 
f|T5|) . {Joir^t + to))to,c can be approximated by 



(j^(r,t + to)X^_^ = j dvUv) ((5(ri(t) -r(t))(5(vi(t) -i;(t))),^_, 

~ j dvUv)f,{r,v,t), (22) 



where fsi^, v, t) is the probabihty density of finding the tracer particle 1 at position r with 
velocity v, which is the same as Eq. fl2Tl) with n{r,t) is replaced by the probability density 



of a tracer particle P{r,t) |l8|. 



From the integration over Vq in Eqs. ( 120]) . we rewrite the current densities as 

{Joir, t + to)>,,^^, ^ u,{r, t + to)P{r, t + to), (23) 

(J^(r,t + to))t,,^^ ^ mn(r,t + to)Mx(?^,t + to)My(^,i + io), (24) 
(JA(r,t + to)>,^^, ^ ^(c/ + 2)n(r,t + to)(r(r,t + to)-T)M,(r',t + to) 

+ ^mn(r, t + to)u^{r, t + to)ucc{r, t + to)- (25) 

Substituting these equations into Eq. (fT9l) . we rewrite the correlation functions in terms of 
the hydrodynamic fields. When the linearized hydrodynamics is used, it is known that cubic 
terms such as the second term in the right hand side of Eq. are negligible [l8j |. 

C. Hydrodynamic equations 

The time evolution of the hydrodynamic fields is described by the following equations 



3^, 



dtn + V ■ (nu) = 0, (26) 

dtu + u- Vu + {nmy^V • n = 0, (27) 

dtT + u-VT + 2{dn)-\U : Vu - AV^T - j^V^n) + TC = 0, (28) 

{dt-DV'')P{r,t) = 0, (29) 

where ( proportional to 1 — is the cooling rate due to inelastic collisions, and D is the 
diffusion constant. We should note that the diffusion equation is introduced to describe the 
diffusion process of the tracer particle 1. 

The pressure tensor U for dilute gases is given by 

Uij = nT5ij - 7] (ViUj + VjUi - ^^ijVkUkj ■ (30) 

Here, the viscosity rj, the heat conductivity A, and the transport coefficient associated 
with the density gradient n, the cooling rate ( and the diffusion coefficient D can be non- 
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dimensionalized as 



where 



V = VoV*, (31) 
A = AoA*, (32) 

fx = n , (33) 

n 

C = — C*, (34) 

^0 

D = ^D*, (35) 
mn 



= ^r(rf/2)7r-V(mT)i/V-('^-i), (36) 
8 



are the viscosity and the heat conductivity in the dilute elastic gas [S^], and i]*, A*, /i*, (*, 
and D* are the constants which depend only on e in dilute cases. 

The hydrodynamic equations (126|1 . (l27j) . (l28!l . and (l29ll have the set of homogeneous 
solutions as 

n(r,t) = n^, (38) 

u{r,t) = 0, (39) 

T(r,t) = Tnit), (40) 

where obeys 

^ = -.H(t)C*TH(t) (41) 

with the characteristic frequency uh = uhTh /riQiTu). Then, we introduce the dimensionless 
deviations of the hydrodynamic variables as 

n{r,t) = nH{l + p{r,t)), (42) 
u{r,t) = UHit)w{r,t), (43) 
T(r,t) = THim+Oir,t)), (44) 

where the variables with the subscript H imply hydrodynamic variables in HCS. It is con- 
venient to introduce the thermal velocity explicitly as 

UH{t) = {TH{t)/mY/\ (45) 



Substituting Eqs. (gSD, (03]), and (gl]) into Eqs. ([26]), ([27]), and ([28]), we obtain a set of 
hydrodynamic equations for the deviation fields p, w, and 9. To avoid to use time dependent 
coefficients, we non-dimensionalize the set of equations by using 



1 



t+to 



dsunis), 



(46) 
(47) 



to 



where Ih = '^UH{t)/vH{t) is the characteristic length in HCS. Here, the dimensionless time 
is proportional to the number of collisions per particle. 

Therefore, the uniform temperature decreases with the dimensionless time as 



TH(t + to) =To(to)exp(-2Cr), 



(48) 



which is nothing but HafF's law 

/ 



351]. We linearize the set of hydrodynamic equations as 



Wk\\ 





—ik 



ik 



9^* d+2 ,.*U2 



a 





—ik 

A* d+2 \ u2 



Wk\\ 



(49) 



r]*k'^ ) Wk_i, 



(50) 

drPk = -k^D*Pk (51) 
in Fourier space, where the Fourier transformed of the hydrodynamic fields are defined by 

fk{r)= [ rf|exp(-zfc-0/(tr). (52) 



Here and w^w are the transversal velocity and the longitudinal component of the velocity, 
which are defined as 



■Wk± = Wk - {k ■ Wk)k, 
Wk\\ = kwk, 

where k = k/\k\. 

The eigenvalues Sa of the matrix in Eq. fl49]) for small k are calculated as 



(53) 
(54) 



d + 2 , „,,3, 



2{d-l 



-X* + 



dC* 



k' + 0{k'), 



~e + oik% 



(55) 



where s+,s_, and s„ correspond to the largest, the smallest and the neutral eigenvalues at 
k = 0, respectively. Here, we note that any sound wave proportional to k does not exist in 
Eq. (ISSj) . This is one of the features in granular gases. In Appendix [Bl we discuss how the 
sound wave disappears in granular gases. In the long-time behavior, we assume that the 
largest eigen-mode dominates the other modes. Thus, the approximate solution of Eq. 
is described as 



Pkit + to) - { + ) exp(s+r) 



Wk\\{t + to) ~ ( — + Wk\\{to) + — lexp(s+r), 

Mt + to) - ^ + J exp(.,r). (56) 

The details of the calculation around here are shown in Appendix \M 
It is easy to find the solutions of Eqs. (!50l) and (ISTIl 

Wk±{t + to) = ^^;fcx(to)e(^*-^'^*'')^ (57) 
Pk{t + to) = Pk{to)e~'"'''\ (58) 

From Eqs. ([38]), (BOD, (05]), (08]), ([MD and ([57D, we obtain the expressions of the 

hydrodynamic fields as 

,^ ( nk{to)k'^ ikuHV^Ukwitp) ^ {d - l)nHTk{U)k'^ \ ^^*_kk^)^ 



. ^ I ik^To{to)nk{to) tkjd - l)Tfc(to) \ ,^2, 

f id-l)Toito)nkito)e , zA:(rf-l)v/mTo(to)^fc||(to) ^ -(C+bfc^ 
°^ - dnnC' dC j"" 

{d - iyTk{to)P ,^,^f^^2)^ 

rf2(*2 ' ^^^^ 

n,x(t + to) = tifc±(to)e-^''*'='^ (62) 

where we introduce Uk = nnPk, Uk = unWk, Tk = TnOk, 

Uk± = Uk- [k ■ Uk)k, (63) 
Uk\\ = kuk, (64) 
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and 

Here, we note that Tk=o decays exponentially as e"''*'^. The origin of the exponential decay 
is the absence of the energy conservation law in granular gases. 

The initial values of the hydrodynamic fields depend on Vq and Tq. In this paper, following 
ref . [l^ , we use the approximate expressions of the initial values of the hydrodynamic fields 

as 

n(r,to) ^ ^^ + ^^(^"-^0(^0)), 
u{r,to) ~ ^W^(r-ro(to)), 

T(r,to) - To(to) + "^"°"^J°^'°V (r-ro(to)), 

P(r,to) ^ Wir-ro{to)). (66) 



D. Results 

Now, let us calculate the correlation functions. Substituting Eqs. (1251) . and fl25]) 
into Eq. f|T9l) with the help of Eqs. fl58|) . fl60|) . fl6Tl) . and flMl) . we can obtain the correlation 
functions. The details of the calculation are shown in Appendices [Cl |Dl and [El 

First, let us calculate (^^(to,^). Substituting Eq. fl23l) into Eq. f|T9l) . we find 

//" dfc 
dvofoito, Vo)vo^ J + to)P-k{t + ^o) 

= C^(to,r) + 4(to,r), (67) 



where C^{tQ,T) and C|)(to,'i") are respectively the contributions from the transverse mode 
and the longitudinal mode for Ukx{t + ^o)- Their explicit definitions are respectively given 
by 

/[ dk 
dvofo{to,vo)vox j j^^Ukx±it + to)P-kit + to), (68) 

r r dk - 

Cj!)(to,r) = Ijf J dvofo{to,vo)vox J j^^Kuk\\{t + to)P-k{t + to). (69) 
The result of CD(tQ,T) (see Appendix lUj) becomes 
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where b is given in Eq. (165|1 . The first term in Eq. (170|) is the contribution from C^{tQ, r), 
and the second term is from C|,(to, t)- We should note that the second term is absent in the 



elastic case 



. The finite contribution from the second term in our problem is related to 



the absence of the sound wave in granular gases. Thus, we obtain the asymptotic behavior 



of C£)ito, t) decaying r '^/^ which is consistent with the simulation by Ahmad and Puri 24 1. 
Second, let us calculate Cr^(tQ,T). Substituting Eq. (l2l|) into Eq. (fT9|) . we rewrite 

Cr,{to,r) ~ —j^ — J dVofo{to,Vo)vo^Voy J J^-^Ukx{t + to)u.ky{t + to) 

= C^^^(to,r) + Cf (to,r) + Cllll(to,r), (71) 



where 



C^^ = m^UHlH / dvQfQ{tQ,Vo)voxVQy / -p---^Ukx±{t + to)u^-ky{t + to), (72) 



dk 

(27) 

/f dk 
dvofoito, Vo)voxVoy J j^^Ukx±it + to)kyU\i^kit + 

dVoMto,Vo)voxVoy J j^-^kxUk\\{t + to)u±_ky{t + tQ), (73) 

/r dk - - 
dvofo{to,Vo)voxVoy J j^-j^kxkyUk\\{t + to)u\\^k{t + to) ■ (74) 

The result of the calculation C^(to,T) becomes (see Appendi:)jDll 

2Toito)\l + a,) 
" d{d + 2)lf, 

""(^^^(4^) +'^(27r(7^* + 26)r) + (s^) )' ^^^^ 

which also obeys C^(to,T) ~ r~^/^. Here, the first term comes from C^-^{tQ,T), the second 
term is from C|^()f:o, t) and the third term is from C|"(to,T). The finite contribution from 
the mixing term C|^(to, t) is also one of the characteristics in granular gases because of the 
absence of the sound wave. 

We should note two important points. First, there should be the logarithmic corrections 
in CnitojT) and Cn(to,T) ior d = 2 through the self-consistent treatment, because the 
coefficients include transport coefficients which show the logarithmic dependence on the 
system size. This is known for the calculation of elastic gases IsgJ. Second, these results 
suggest that the diffusion constant and the viscosity may depend on the system size in 
two-dimensional systems. We discuss these system size dependences briefly in Sec. HVl 
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In final, let us evaluate Cx(tQ,T). Substituting Eq. (125|) into Eq. (fT9|) . we obtain (see 
Appendi jEj) 

{d + 2)nH f,ru ^ /"^^o (rf + 2)To(to) 



<^A(to,i-) ^ — y o?Vo/o(io,t'o)fox 1^- 



2 



X J -^^Uk^it + to)T^k{t + to) 

Ci{to,T) + Cl{to,T), (76) 



where 



<^a(^o,^) = / fivo/o(to,Wo)wox ( ^ ^ 

X J j^Uk^±{t + to)T.kit + to), (77) 

C\(^o,r) = ^ / dvofo{to,Vo)vo^ {— 1 

/dk - 
j^Kukii{t + to)T_k{t + to). (78) 

The calculation of Cx{to, r) gives 

7r(^ + 2)^(rf-l)ro(to)M {d-l)e-i'^ , 2de-C'- ^ 
2d^mC*n% y (27r(r/* + 26)r)— (87r6r) — / 

where we introduce A as 

A = 2(rf- 1) + af(M- 10). (80) 

Here the first term in the right hand of Eq. fl79|) comes from Cj;{tQ, r), and the second term 
in Eq. (ff9l) is from C]I(to, t). The result in Eq. ( ff9l) is characteristic one for granular fluids, 
because CA(to)''") has the fast decay obeying r^^'^+^^/^e"^*'^. The absence of the long-time 
tail in C\(tQ, r) is the result of the absence of the energy conservation law. 

In our calculation, we assume the linearity of the fluctuations around the homogeneous 
cooling state. This assumption might be invalid to describe the long-time behavior of freely 
cooling granular gases. Indeed, it is well-known that the homogeneous cooling state is 
unstable, and the system develops into inhomogeneous complicate patterns. Hence, it is not 
clear whether our analytic calculation of the correlation functions is valid in the wide range 
of time evolution of freely cooling states. 
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Before we close this section, let us comment on the initial condition dependence of our 
result. Through our calculation, we assume that the average is taken over possible configu- 
ration at time to- Although we introduce /; in Eq. (pTi) . we use it only to derive Eqs. 
(^^, and (^^. Fortunately, our analysis for Cr)itQ,t), C^(to,t), and Cx{tQ,t) suggests that 
all correlation functions can be factorized. Namely, the correlations are represented by the 
products of functions of to and functions of t. It is obvious that the important point for long 
time behavior is the functions of t, and long-time tails are function of t. 



III. NUMERICAL SIMULATIONS 



To verify the validity of our theory, in this section, we perform two-dimensional simula- 
tions to calculate Coito, r) and CA(to, t), where we only calculate the kinetic part of CA(to, r) 
in Eqs. and (jTj). Here, we do not show any result of C^(to,T) because there are large 
fluctuations in the numerical data and the results are not characteristic one. To simulate 
the system, we employ the event driven method for hard-core particles in which the collision 
rule is given by Eq. ([1]), where the scaled time by the collision frequency is the only relevant 
time. The actual code of the simulations is based on the efficient algorithm developed by 
Isobe 4J]. We set the volume fraction u = z/c/2, where Uf. is the close packing fraction of 



hard disks. We choose the equilibrium state as the initial state at t = 0. In addition, we fix 
the time when we start the measurement as zero. Therefore, in this section, we replace the 
notation Ca(to,r) by Cair). In our simulation, To(to), m and a are set to be unity, and all 
quantities are converted to dimensionless forms. We should note that our numerical system 
is not dilute which is in contrast to the assumption of the theory. However, as mentioned in 
the previous section, we expect that the contribution from the potential term is not essen- 
tially important when the density is enough lower than the jamming transition point. Thus, 
we regard the average volume fraction u = z/c/2 is enough "dilute". 

In our simulation, the number of the particles we use is 65536 for the calculation of 
with the ensembles of 100 different initial conditions for the case e = 1.0 and 10 initial 
conditions for the case e < 1.0. On the other hand, to suppress the large fluctuation in the 
calculation of Cx{t), we use a smaller system with = 1024 with the ensemble of 2250000 
initial conditions for the case e = 1.0 and 400000 initial conditions for the case e = 0.95 {s?]. 
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A. The velocity auto-correlation function 



Figure [U shows the behavior of Cd{j) as the function of r with e = 1.0,0.9, and 0.8. 
For e = 1.0, C£,{t) approximately obeys r^^ as indicated in the previous studies |l^. For 
granular gases with e < 1.0, Cj:,{t) decays as when the time r is smaller than a certain 
threshold value Tp. This behavior in the early regime is consistent with the result of our 
theory. On the other hand, when r > rp, CdIt) seems to decay faster than r^^. This 
behavior is contrast to the prediction of our theory. It is not surprising that our theory 
cannot be used for r > rp. This is because the system is no longer in HCS, and we 
conjecture that our theory is valid in HCS. The fact that the system is not in HCS for 
r > Tp is confirmed from Fig. [2|^a). This figure shows the snapshot of the system with 
e = 0.8 around r = 22 which is nearly equal to rp. From Fig. [2t^a), we find the existence of 
inhomogeneous clusters. We also note that the threshold value rp decreases as e decreases. 




FIG. 1: (a) Co{t) as the function of r with e = 1.0,0.9. (b) Cd{t) as the function of r with 
e = 0.8. th and tl are the time of the violation of Haff's law and the time starting to obey 
T(r) ~ r~^, respectively. 



We also calculate T(r) as the function of r for e = 0.8 (Fig. [3]). We confirm that the 



temperature decreases in terms of Haff's law fl48p 35l | for r < 22. However, when r > 22, 
Haff's law is no longer valid and the system becomes inhomogeneous. Here, we introduce 
Th as the time when Haff's law is violated. Our picture that the violation of the theoretical 
description is related to the formation of inhomogeneous clusters can be verified through the 
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(a) 



(b) 



FIG. 2: (a) The snapshot of the density field at r = 22. (b) The snapshot of the density field at 
r = 70. Both figures are obtained from simulation of e = 0.8. 

numerical comparison of Tp with th (FigJlj), in which th is close to rp for all e. This fact 
supports that our theory is valid in HCS, but the theory may not be used when the system 
goes into the inhomogeneous state. 

T 

0.1 
0.01 

0.001 
0.0001 

1 10 100 1000 

T 

FIG. 3: T{t) as the function of r with e = 0.8. See the caption of Fig. [T]for the definitions of th 
and Tl- 

From the data shown in Fig. [H we find another interesting behavior. When the system 
goes into the late regime, Cd{t) seems to recover the theoretical behavior: Cd{t) ~ t~^. 
Indeed, the simulation for e = 0.8 clearly shows the existence of two regions obeying as 
shown in Fig. [T](b). The time when Cz)(r) recovers the theoretical behavior decreases as e 
decreases. 
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FIG. 4: th and rp as the function of 1 — 



It is interesting that the numerical result becomes consistent with the theoretical predic- 
tion in the late regime. We do not have any definite answer of this mysterious behavior. 
However, we note that the linearized fiuctuation hydrodynamics is valid even in the late 
stage 12|, and the decay of temperature T{t) ~ t"''/^ in the late stage can be calculated 
from the fiuctuation of the linearized hydrodynamics [38] . Here, we have confirmed that the 
scaled r is almost proportional to the real time in the late stage, which is contrast to the 
relation in HCS. Let us introduce ti as the time when T(r) starts to obey r~'^/^, as shown in 
Fig. [3l From the comparison of and the behavior of CdIt) in Fig. [T](b), we find the time 
that Cd{t) starts to recover the theoretical behavior is nearly equal to r^. We also suppose 
that the theoretical behavior recovers because the density becomes uniform in clusters in 
the late stage. Indeed, the uniform density inside clusters may be verified in Fig. [2](b) at 
r = 70 which is nearly equal to r^. 



B. The correlation function for the heat flux 



Figure [5] shows the result of Cx{t) as the function of r for e = 1.0 and 0.95. From this fig- 
ure, we find that Cx{t) decays with an exponential function of r when e = 0.95, while Cx{t) 
has a long tail obeying Cx ~ r^^ for e = 1.0 though the range to obey is short. This 
behavior is different from the behavior of Cd{t) as expected from the theoretical prediction. 
The fluctuation of the data for Ca(t), however, is too large to verify the quantitative validity 
of the theoretical prediction Ca(t) ~ _ T~('^+2)/2e~^*'^. It should be noted that our theory 
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suggests negative C\{t) but the simulation does not show negative Ca(t) 
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Cx(t) ~ exp(-0.07*T)~^. ^^"^ ~~ - 
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FIG. 5: Cx{t) as the function of r with e = 1.0, 0.95. 



IV. DISCUSSION AND CONCLUSION 

In this paper, we analytically calculate the behavior of the time correlation functions 
CD(to,T), Cr^(to,T), and Cx{to,T). From our calculation, we predict that CoitoyT) and 
Crj{to,T) obey r"*^/^, while Cx{tQ,T) decays as r^(^+'^y^e~''*'^ in the freely cooling system. 
Although the behavior of Cx{t(),T) is different from that in elastic gases, the behaviors 
of C£){tQ,T) and C^(to,T) are similar to those in the elastic cases. We also perform the 
numerical simulation of the freely cooling system to verify the results of the theoretical 
prediction. Through our numerical calculation of correlation functions, we find that our 
prediction is valid in HCS. Although our prediction cannot be used in the middle regime of 
inhomogeneous state, the theoretical behavior is recovered in the late regime. 

The prediction of our theory is completely different from Kumaran's prediction I23!] for 
the sheared granular fluids in which the correlation functions obey t~^'^^'^. One reason of the 
difference between the theories comes from the difference of the base states in the analyses. 
In our theory, we choose the homogeneous cooling state as the base state, while the base 
state in the sheared granular flow is the uniform shear flow in which the velocity gradient is 
uniform, the temperature and the density are uniform. 

In the elastic gases, the autocorrelation functions are directly related to the transport 
coefficients. Actually the diffusion coefficient, the viscosity and the heat conductivity are 
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respectively given by D oc dtCnit), f] oc dtC^{t), and A oc dtC\{t) for conventional 
: luids. It is known that such simple Green-Kubo formula may not be valid in granular gases 
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43|. 



The anomalous relation between the transport coefficient and the correlation functions in 
granular fluids is easily verified from the following simple argument. If we assume Ci){to,T) = 
{v(t + to) ■ '^(to))to ~ ^o(^o)0(t), then, we obtain the relation 

for sufficiently large r, where (r^(r)) is the mean square displacement of the tracer particle 
during r. Although the explicit expression seems to be different from the conventional 
Green-Kubo formula, the logarithmic singularities in two-dimensional systems are still valid 
at least near e = 1. In fact, the integral can be evaluated as 0(r ds(j){s)) = 0(r log r) for 
d = 2 and r < l/{2(*) with l/{2(*) is proportional to (1 — e^)~^. For this case we might 
replace the cut-off time by the characteristic time depending on the linear system size L such 
as L'^/D, and D may be determined by the self-consistent way. Thus, we expect that the 
diffusion coefficient depends on the system size as the result of the logarithmic singularity. 
We also indicate that Eq. flSTl) implies that the diffusion disappears in the limit of large r, 
because the exponential term is negligible. This is also consistent with the observation in 
our simulation. 

In this paper, we adopt Eq. f|T7|) as the initial weight function which depends on the 
initial temperature. However, this dependence may be removed if we map HCS onto a 
steady state 12]. We expect that this mapping make our calculation clearer. 

We also assume that the coefficient of restitution e is a constant in contrast to actual cases 
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261,121 



particles 



It may be possible to extend our results to the gases of viscoelastic 
47( 1 . However, we believe that our prediction for the algebraic time dependence 



of the correlation functions is unchanged. 



3^. The 



We adopt, here, the phenomeno logical approach developed by Ernst et al. 
advantages of this method is that the calculation is simple and the physical pictures i.e. 
the roles of conserved quantities are clear. This approach is successful for the first step to 
indicate the relevant roles of long-time tails. However, the approach has the limitation in 
which our theoretical prediction cannot be used in the middle stage. To improve this, we 



may need a more fundamental approach starting from the Liouville equation 
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48| 



In conclusion, we confirm the existence of tlie long-time tail in the auto-correlation func- 
tions for the velocity and the shear stress, while it does not exist for the heat flux. These 
results are consistently obtained from the theory and simulation. Thus, the transport coef- 
ficient based on inelastic Enskog equation should be modified for two-dimensional systems. 

One of the authors (HH) thanks J. Wakou for his useful comments. We also appreciates 
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Japan (Grant No. 18540371) and the Grant-in- Aid for the 21st century COE "Center for 
Diversity and Universality in Physics" from MEXT, Japan. The numerical calculations were 
carried out on Altix3700 BX2 at YITP in Kyoto University. 



APPENDIX A: THE TIME EVOLUTION OF THE HYDRODYNAMIC FIELDS 



In this appendix, we calculate the solution of Eq. (H^ . The eigenvalue s of the matrix 
in Eq. fH9!) satisfies the equation 



-s[s-C + ^V*k" ][s + C + 



d + 2 



-S-C + 



-A A;^ 



2{d-l] 



0. 



d ^ ^ ' 2(d- 1) 2{d-l) 
To discuss hydrodynamic behaviors, we expand the eigenvalue around /c = as 

S = So + Sik + 52^^ H . 



(Al) 



(A2) 



Substituting this equation into Eq. (lAip . we find that there exists three eigenvalues s+, s_, 
and Sn as Eqs. fl55|) . Introducing the eigenvector corresponding to the eigenvalue 
with the condition \Xa\'^ = 1, we find that the solution of Eq. (H9l) is expressed as 



^ pk{t+to) ^ 

Wk\\{t + to) 
\ dkit + to) J 

where we define 



a+ exp(s+r)X+ + a_ exp(s_r)X_ + a„ exp(s„r)X„, (A3) 



(A4) 
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Near k = 0, the eigen- vector X+ can be represented by 



X 



{d-l)ik 



l + Oik) 



(A5) 



and 



Pkito)k \ l)ik9k{to) 

+ Wkiato) H — H 



(A6) 



(* ' ""'^"^""^ ■ dC* 

When the time r is enough large, the mode corresponding to the largest eigenvalue s+ is 
only the relevant mode. Hence, we obtain the approximate solution of Eq. ( l49l) as Eq. ( 156|1 
in the long-time behavior. 



APPENDIX B: THE ABSENCE OF THE SOUND WAVE 

It is well known that the eigenvalue s has the term proportional to k which corresponds 
to the sound wave for elastic gases. However, as shown in Eq. fl55|) . for granular gases with 
e < 1.0, the eigenvalue s does not have the term proportional to k. This might be puzzling 
because the term proportional to k does not appear in the elastic limit of Eq. ( !55l) at the 
first glance. 

In order to understand this puzzle, let us show the relevant point of the calculation of s 
in Appendix |X1 Substituting Eq. flA2l) into flAl|) with equating coefficients of each power of 
/c, we obtain the expression of s. 

In the 0-th and first order of /c, we obtain the equation 

So(so- 0(^0 + = 0, (Bl) 

s\{^sl-C')=Q. (B2) 

From Eq. (IBip . we obtain 

So = aC\ (B3) 

where a = 0, ±1 with oc 1 — e^. Substituting this equation into Eq. ( ]B2[) . we find 

s?(3a^ - 1)C*^ = 0. (B4) 
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For (* 7^ 0, thus, Si only satisfies 

si = 0. (B5) 

However, for tlie elastic case with (* = 0, si is not determined by this equation but by the 
equation obtained from the coefficients of in Eq. flAip . and we find that si has some 



non-zero values. Thus, the elastic limiti is rather singular in granular fluids. This is the 
reason that the there are no sound wave when e < 1.0. 

APPENDIX C: THE CALCULATION OF CdC^o t) 

In this appendix, we calculate the correlation function Eq. (170]) . From Eq. (1^ . Cr){tQ, r) 
is expressed as 

Coito, r) = C^to, r) + 4 (to, r), (CI) 

where C^{to,T) and C^j){to,T) are defined in Eqs. ( l68l) and ( l69l) . 

Substituting Eqs. fl58l) and flM|) into Eq. fl68|) . we obtain the expression of C^(tQ,T) as 

Ch = j dvofo{to,vo)vo. I ^«fc,a(to)P-.(to)e-(^''*+^*)'='^ 



uhIh 



(d-l)To(to)/ 1 



dmljj \2Tc{ri* + 2D*)t ^ 
where we use the relation derived from Eq f|T3|) as 



(C2) 



lim\Wk\ = l. (C3) 



In addition, we also use the relation 

J dkkik = k/d- (C4) 
Similarly, substituting Eqs. fl58l) and fl60l) into Eq. fl68|) with the use of Eq. flC3|) . we 
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obtain the expression for Cj^ito, r) as 



To (to) f dk 



mljj J {27iY " 
Toito) f 1 



dmlfj \A7r (b + D*) 
From Eqs. ^ and we obtain Eq. ([70 



(C5) 



APPENDIX D: THE CALCULATION OF Cr,{to,T) 

Let us calculate C^(to,T). From Eq. (ITT]) . C^(to,T) is expressed as 

C,(to, r) = C^^^(to, r) + (to, r) + ClJ" (to, r), (Dl) 

where C;j^-'-(to, t), C,^"(to,''") and C|"(to,T) are respectively defined by Eqs. (1721) . (173|l and 
dZH). 

Substituting Eq. (1641) into Eq. (1721) . we obtain the expression of C^-'-(to,T) as 

f f dk 

C^^(to,r) = w?nHljf I dvofo{to,vo)vo^voy / 77;;;;^Ufcx±(io)^i±-fes,(^o)e"''*'''^ 



J dVofo{to,Vo)VoxVoy J j^^ivox - K^Q ' k){voy - kyVo ■ k)\Wk\ 

(<i2 - i)r„((„) 



where we use Eq. ( IC4p and 



CLKKiKjKiKm — _|_ 2) ' V 



Similarly, substituting Eqs. (1601) and (IMI) into Eq. (173|) . we obtain the expression of 
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C|^(to,T) as 

/f dk 



/f dk 
dVofo{to,Vo)VoxVOy J Jj^kxUk\\{to)Uk±y{to) 



dvofo{to,vo)voxVoy / jTT^ivox - K{vo ■ k))ky{vo ■ k)\Wk\'^e'^^'^'''^''r ^ 



uhI'hJ V (2vr) 



TTt^ Tx f (* dk 

^mnn y dvofo{to,vo)voxVoy J j^kivo ■ k){voy - kyivo ■ k))\Wk\^e-(^'''+'y^ 



I'h J (2vr) 



d 



'^°^^°)\l + af)f--i--) . (D4) 



{d + 2)ljj^ ^ ^ \2tt{7]* + 2b)T 
Substituting Eq. (pUj) into Eq. (171|1 . we obtain the expression of C|"(ito,T) as 

dvofo{to, Vo)voxVoy J j^--^kxkyUk\\{tQ)u\\^k{tQ) 



—jj- J dvofo{to,vo)voxVoy J j^^kky{vo- kflWkl'^e'^''''^'' 

7o(to)^(l + a^) f dk opplry |2^-2bfcV 



From these equations, we obtain Eq. ([75 



APPENDIX E: THE CALCULATION OF Cx{to,T) 

In this Appendix, we calculate C\(to,T). From Eq. (176|) . CA(to,''") is the combination of 
two modes. 

Cxito, r) = Ciito, r) + Clito, r), (El) 

where C^{to,T) and c|(to,T) are respectively defined by Eqs. fl77j) and flTHj) . Since the 
calculation is complicated, we separate the contribution from the transverse mode and that 
from the longitudinal mode in the explanation. 
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1. The contribution from the transverse mode 



Let us evaluate C^{to,T) which is the contribution from the transverse mode at first. 
Substituting Eqs. flM]) and flM|) into Eq. fl77|) . we obtain the expression of C^{to,T) as 

\d - l)To{t^)nu{ta)e ^ ik{d - l)v/mTo(to)M,.||(to) {d - lfTu{t^)k^ 



dnnC*^ dC* d'^C*^ 

Ci^{to,T) + Ci^{t,,T), (E2) 



where we introduce 



(27r) 



7^u..^(to)n_.(to)A:^e-{^+(^''*-^^)'=^}^ (E4) 

Here, we use the fact that the term proportional to Ukx±{tQ)u-k\\{tQ) in Eq. (161]) does not 
survive in the integration over vq. 

Substituting Eq. (!66|) into Eq. (lE3p . we obtain 



+ 2)^(. - l)%(to)3 ^ ^^^^ ^ 10)) _ --^1 (E5) 



2d^mC% {2TT{r]* + 2b)Ty 
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Substituting Eq. (166|) into Eq. (lE4p . we obtain 



Ti{d + 2)2(rf - l)2To(to)3af / e-^*^ 



Substituting Eqs. flESp and flE6p into Eq. flE2p . we obtain 



(E6) 



where A is defined by Eq. fl80p . 

2. The contribution from the longitudinal mode 

Substituting Eqs. flBUl) and flM]) into Eq. fl75]) . we obtain the expression of C|(to,T) as 

C^!|(to,r) = +1^^ I dWo(to,t;o)t^o.Ho-(rf + 2)To(to)) / (|^^xe-{^+^^^^}- 

/ ik^/%(to)nk{to) , , , l)v^Tfc(to) 
X t;^ ^ + ^ifcil(^o) + 



{d - l)To(to)nfc(to)P zA;(rf-l)v/mTo(to)nfc||(to) _ (rf - l)^Tfc(to)A;^ 
Cf (to,r) + Cf(to,r), (E8) 



where we introduce 



(to,r) = JA±3^—^^ I dvofo{to,vo)vo4mv',-{d + 2)To{to)) 



2d^C*% 

X J ^i«fc||.(to)T_.(to)^'e-{^+^^'='}^ (E9) 



dk 



^^^yKu^ito)n.kito)k'e~i'^^'''"}\ (ElO) 
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Substituting Eq. fl66|l into Eq. flE9|) . we obtain the expression for Cf'^{to,T 



IIT,, (rf + 2)(rf-l)2 



Mvo ■ k){mvl - dTo(to))|W^fc|'A;V{^+(5''*+'')'='h 



dk 



2d^mC*% J {2nY " 

nid + 2)\d - ifTM^ ,^ „ . . , . / e-^*- 



-(2 + a^(rf+10)) — --^ . (Ell) 



For C|"(to,T), by substituting Eq. (166!) into Eq. (lElOp . we obtain 

'-^A (^o,r) = — — / dvQto[h,VQ)vQa^[mvQ- [d + 2)TQ[tQ)) 



2dnHC*% 

{d + 2f{d-l)T,{t,fa^ f J^p, 2 ^{(+2bk^}. 
2dmC*% J (27r)<^ " 

7r(rf + 2)2(rf-l)2To(to)'af / e'^* 



ii+2 



Thus, from Eqs. flETTl) and flET2l) we obtain C|(to,r) as 



(E12) 



rWu ^^ 7r(rf + 2)^(rf-l)^To(to)M / e'^'- \ 
Finally, from Eqs. flE7p and flElSp associated with Eq. (1751) . we obtain Eq. (jTHD- 
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